Data analysis for:
Padilla Perez et al., (2022). The correlated evolution of foraging mode and reproductive output in lizards. Proceedings of the Royal Society B. Dylan J. Padilla Perez, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.Library
require(AICcmodavg)
require(ape)
require(caper)
require(geiger)
require(kableExtra)
require(MuMIn)
require(nlme)
require(phytools)
require(plotrix)
require(shape)
require(car)R.version _
platform x86_64-pc-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 4
minor 1.0
year 2021
month 05
day 18
svn rev 80317
language R
version.string R version 4.1.0 (2021-05-18)
nickname Camp PontanezensessionInfo() R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS: /packages/7x/R/R-4.1.0/lib/libRblas.so
LAPACK: /packages/7x/R/R-4.1.0/lib/libRlapack.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] car_3.0-13 carData_3.0-5 shape_1.4.6 plotrix_3.8-2
[5] phytools_1.0-1 maps_3.4.0 nlme_3.1-153 MuMIn_1.43.17
[9] kableExtra_1.3.4 geiger_2.0.7 caper_1.0.1 mvtnorm_1.1-3
[13] MASS_7.3-55 ape_5.6-2 AICcmodavg_2.3-1 knitr_1.37
[17] rmarkdown_2.13
loaded via a namespace (and not attached):
[1] Rcpp_1.0.8.3 subplex_1.7 svglite_2.1.0
[4] lattice_0.20-45 digest_0.6.29 unmarked_1.1.1
[7] R6_2.5.1 plyr_1.8.7 stats4_4.1.0
[10] evaluate_0.15 coda_0.19-4 httr_1.4.2
[13] rlang_0.4.12 rstudioapi_0.13 raster_3.5-15
[16] jquerylib_0.1.4 phangorn_2.8.1 Matrix_1.4-0
[19] combinat_0.0-8 splines_4.1.0 webshot_0.5.3
[22] stringr_1.4.0 munsell_0.5.0 igraph_1.2.11
[25] compiler_4.1.0 numDeriv_2016.8-1.1 xfun_0.30
[28] systemfonts_1.0.4 pkgconfig_2.0.3 mnormt_2.0.2
[31] tmvnsim_1.0-2 htmltools_0.5.2 expm_0.999-6
[34] quadprog_1.5-8 codetools_0.2-18 viridisLite_0.4.0
[37] grid_4.1.0 lifecycle_1.0.1 jsonlite_1.8.0
[40] xtable_1.8-4 magrittr_2.0.3 scales_1.1.1
[43] stringi_1.7.6 scatterplot3d_0.3-41 sp_1.4-6
[46] xml2_1.3.3 bslib_0.3.1 fastmatch_1.1-3
[49] deSolve_1.30 tools_4.1.0 glue_1.6.2
[52] abind_1.4-5 parallel_4.1.0 fastmap_1.1.0
[55] survival_3.3-1 yaml_2.3.5 colorspace_2.0-3
[58] terra_1.5-21 rvest_1.0.2 VGAM_1.1-5
[61] clusterGeneration_1.3.7 sass_0.4.1evo_bd <- read.csv("evo_body_size_mass2.csv", row.names = 1)
evo_bd <- evo_bd[evo_bd$diet != "Herbivorous", ]
evo_bd <- evo_bd[evo_bd$Family != "Sphenodontidae", ]
tree <- read.nexus("time.tree.nex")
str(evo_bd) 'data.frame': 588 obs. of 25 variables:
$ reproductive.mode : chr "Oviparous" "Oviparous" "Oviparous" "Oviparous" ...
$ Genus : chr "Acanthocercus" "Acanthodactylus" "Acanthodactylus" "Acanthodactylus" ...
$ Family : chr "Agamidae" "Lacertidae" "Lacertidae" "Lacertidae" ...
$ female.SVL : num 113.9 64.9 60.1 68.6 61 ...
$ hatchling.neonate.SVL : num 34 25 26.5 29 30 35.8 24.3 65.3 47 27 ...
$ female.log : num 2.1 1.8 1.8 1.8 1.8 1.8 1.7 2.1 2 1.9 ...
$ hatch.log : num 1.5 1.4 1.4 1.5 1.5 1.6 1.4 1.8 1.7 1.4 ...
$ intercept : num -4.69 -4.54 -4.54 -4.54 -4.54 ...
$ slope : num 3.1 2.95 2.95 2.95 2.95 ...
$ female_mass : num 349 187 173 198 176 ...
$ hatchling_mass : num 100.9 69.2 73.7 81 84 ...
$ femass.log : num 1.7 0.8 0.7 0.9 0.7 0.9 0.5 1.6 1.4 1.3 ...
$ hatcmass.log : num 0.07 -0.42 -0.34 -0.23 -0.18 0.04 -0.45 0.74 0.24 -0.24 ...
$ rawFem : num 50.06 6.38 5.09 7.52 5.32 ...
$ rawHatc : num 1.17 0.38 0.45 0.59 0.65 1.1 0.35 5.44 1.74 0.57 ...
$ largest.clutch : int 17 8 8 8 7 8 6 1 2 18 ...
$ smallest.clutch : int 1 3 1 1 2 1 1 1 1 8 ...
$ largest.mean.clutch.size : num 11.4 5 6 4.4 4.8 3 3.3 1 2 14 ...
$ smallest.mean.clutch.size : num 11.4 4.8 2.6 3 2.5 2.6 2.5 1 2 11.5 ...
$ Longitude.centroid..from.Roll.et.al..2017.: num 30.3 34.8 19.8 -3.2 25.1 ...
$ Latitude.centroid..from.Roll.et.al..2017. : num -9.53 31.17 23.35 37.24 31.05 ...
$ Activity.time : chr "Diurnal" "Diurnal" "Diurnal" "Diurnal" ...
$ foraging.mode : chr "sit and wait" "widely foraging" "widely foraging" "mixed" ...
$ diet : chr "Carnivorous" "Carnivorous" "Carnivorous" "Carnivorous" ...
$ Binomial : chr "Acanthocercus_atricollis" "Acanthodactylus_beershebensis" "Acanthodactylus_boskianus" "Acanthodactylus_erythrurus" ...## Scaled mass index for females
par(mar = c(6, 6, 1, 1), mgp = c(4, 1, 0))
plot(rawFem ~ female.SVL, data = evo_bd, ylab = "Female mass (g)", xlab = "Female length (mm)", las = 1)
par(mar = c(5.1, 4.1, 4.1, 2.1), mgp = c(3, 1, 0))
plot(log(rawFem) ~ log(female.SVL), data = evo_bd, ylab = "Female mass ln(g)", xlab = "Female length ln(mm)", las = 1)
step1 <- lm(log(rawFem) ~ log(female.SVL), data = evo_bd)
summary(step1)
Call:
lm(formula = log(rawFem) ~ log(female.SVL), data = evo_bd)
Residuals:
Min 1Q Median 3Q Max
-2.6266 -0.1464 -0.0359 0.2153 0.7334
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.51927 0.11602 -90.67 <2e-16 ***
log(female.SVL) 2.96706 0.02686 110.48 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3334 on 586 degrees of freedom
Multiple R-squared: 0.9542, Adjusted R-squared: 0.9541
F-statistic: 1.221e+04 on 1 and 586 DF, p-value: < 2.2e-16L0 <- mean(evo_bd$female.SVL)
Mi <- evo_bd$rawFem
Li <- evo_bd$female.SVL
evo_bd$M_females <- Mi*(L0/Li)^3
## Scaled mass index for hatchlings
plot(rawHatc ~ hatchling.neonate.SVL, data = evo_bd, ylab = "Hatchling mass (g)", xlab = "Hatchling length (mm)", las = 1)
plot(log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd, ylab = "Hatchling mass log(g)", xlab = "Hatchling length log(mm)", las = 1)
step1h <- lm(log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd)
summary(step1h)
Call:
lm(formula = log(rawHatc) ~ log(hatchling.neonate.SVL), data = evo_bd)
Residuals:
Min 1Q Median 3Q Max
-1.70359 -0.30272 0.02864 0.17771 2.09892
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -10.11741 0.14347 -70.52 <2e-16 ***
log(hatchling.neonate.SVL) 2.82962 0.04232 66.87 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.4173 on 586 degrees of freedom
Multiple R-squared: 0.8841, Adjusted R-squared: 0.8839
F-statistic: 4471 on 1 and 586 DF, p-value: < 2.2e-16L0h <- mean(evo_bd$hatchling.neonate.SVL)
Mih <- evo_bd$rawHatc
Lih <- evo_bd$hatchling.neonate.SVL
evo_bd$M_hatchlings <- Mih*(L0h/Lih)^2.8
head(evo_bd) reproductive.mode Genus Family
Acanthocercus_atricollis Oviparous Acanthocercus Agamidae
Acanthodactylus_beershebensis Oviparous Acanthodactylus Lacertidae
Acanthodactylus_boskianus Oviparous Acanthodactylus Lacertidae
Acanthodactylus_erythrurus Oviparous Acanthodactylus Lacertidae
Acanthodactylus_pardalis Oviparous Acanthodactylus Lacertidae
Acanthodactylus_schreiberi Oviparous Acanthodactylus Lacertidae
female.SVL hatchling.neonate.SVL female.log
Acanthocercus_atricollis 113.9 34.0 2.1
Acanthodactylus_beershebensis 64.9 25.0 1.8
Acanthodactylus_boskianus 60.1 26.5 1.8
Acanthodactylus_erythrurus 68.6 29.0 1.8
Acanthodactylus_pardalis 61.0 30.0 1.8
Acanthodactylus_schreiberi 69.4 35.8 1.8
hatch.log intercept slope female_mass
Acanthocercus_atricollis 1.5 -4.686 3.105 349.0
Acanthodactylus_beershebensis 1.4 -4.543 2.951 187.0
Acanthodactylus_boskianus 1.4 -4.543 2.951 172.8
Acanthodactylus_erythrurus 1.5 -4.543 2.951 197.9
Acanthodactylus_pardalis 1.5 -4.543 2.951 175.5
Acanthodactylus_schreiberi 1.6 -4.543 2.951 200.3
hatchling_mass femass.log hatcmass.log rawFem
Acanthocercus_atricollis 100.9 1.7 0.07 50.06
Acanthodactylus_beershebensis 69.2 0.8 -0.42 6.38
Acanthodactylus_boskianus 73.7 0.7 -0.34 5.09
Acanthodactylus_erythrurus 81.0 0.9 -0.23 7.52
Acanthodactylus_pardalis 84.0 0.7 -0.18 5.32
Acanthodactylus_schreiberi 101.1 0.9 0.04 7.78
rawHatc largest.clutch smallest.clutch
Acanthocercus_atricollis 1.17 17 1
Acanthodactylus_beershebensis 0.38 8 3
Acanthodactylus_boskianus 0.45 8 1
Acanthodactylus_erythrurus 0.59 8 1
Acanthodactylus_pardalis 0.65 7 2
Acanthodactylus_schreiberi 1.10 8 1
largest.mean.clutch.size
Acanthocercus_atricollis 11.4
Acanthodactylus_beershebensis 5.0
Acanthodactylus_boskianus 6.0
Acanthodactylus_erythrurus 4.4
Acanthodactylus_pardalis 4.8
Acanthodactylus_schreiberi 3.0
smallest.mean.clutch.size
Acanthocercus_atricollis 11.4
Acanthodactylus_beershebensis 4.8
Acanthodactylus_boskianus 2.6
Acanthodactylus_erythrurus 3.0
Acanthodactylus_pardalis 2.5
Acanthodactylus_schreiberi 2.6
Longitude.centroid..from.Roll.et.al..2017.
Acanthocercus_atricollis 30.28
Acanthodactylus_beershebensis 34.83
Acanthodactylus_boskianus 19.84
Acanthodactylus_erythrurus -3.20
Acanthodactylus_pardalis 25.10
Acanthodactylus_schreiberi 33.23
Latitude.centroid..from.Roll.et.al..2017.
Acanthocercus_atricollis -9.53
Acanthodactylus_beershebensis 31.17
Acanthodactylus_boskianus 23.35
Acanthodactylus_erythrurus 37.24
Acanthodactylus_pardalis 31.05
Acanthodactylus_schreiberi 35.05
Activity.time foraging.mode diet
Acanthocercus_atricollis Diurnal sit and wait Carnivorous
Acanthodactylus_beershebensis Diurnal widely foraging Carnivorous
Acanthodactylus_boskianus Diurnal widely foraging Carnivorous
Acanthodactylus_erythrurus Diurnal mixed Carnivorous
Acanthodactylus_pardalis Diurnal sit and wait Carnivorous
Acanthodactylus_schreiberi Diurnal widely foraging Carnivorous
Binomial M_females
Acanthocercus_atricollis Acanthocercus_atricollis 22.23411
Acanthodactylus_beershebensis Acanthodactylus_beershebensis 15.31746
Acanthodactylus_boskianus Acanthodactylus_boskianus 15.38844
Acanthodactylus_erythrurus Acanthodactylus_erythrurus 15.28782
Acanthodactylus_pardalis Acanthodactylus_pardalis 15.38233
Acanthodactylus_schreiberi Acanthodactylus_schreiberi 15.27571
M_hatchlings
Acanthocercus_atricollis 0.9884380
Acanthodactylus_beershebensis 0.7593747
Acanthodactylus_boskianus 0.7638861
Acanthodactylus_erythrurus 0.7781121
Acanthodactylus_pardalis 0.7796118
Acanthodactylus_schreiberi 0.8043124check <- name.check(tree, evo_bd)
rm_phy <- check$tree_not_data
rm_dat <- check$data_not_tree
ctree <- drop.tip(tree, rm_phy)
cdat <- subset(evo_bd, subset = evo_bd$Binomial %in% ctree$tip, select = names(evo_bd)[1:27])
name.check(ctree, cdat) [1] "OK"## Mean clutch size
cdat$mclutch <- (cdat$smallest.mean.clutch.size + cdat$largest.mean.clutch.size) / 2
## Reproductive output
cdat$rep_out <- (cdat$mclutch * cdat$M_hatchlings)
str(cdat) 'data.frame': 485 obs. of 29 variables:
$ reproductive.mode : chr "Oviparous" "Oviparous" "Oviparous" "Oviparous" ...
$ Genus : chr "Acanthocercus" "Acanthodactylus" "Acanthodactylus" "Acanthodactylus" ...
$ Family : chr "Agamidae" "Lacertidae" "Lacertidae" "Lacertidae" ...
$ female.SVL : num 113.9 64.9 60.1 68.6 61 ...
$ hatchling.neonate.SVL : num 34 25 26.5 29 30 35.8 24.3 65.3 47 27 ...
$ female.log : num 2.1 1.8 1.8 1.8 1.8 1.8 1.7 2.1 2 1.9 ...
$ hatch.log : num 1.5 1.4 1.4 1.5 1.5 1.6 1.4 1.8 1.7 1.4 ...
$ intercept : num -4.69 -4.54 -4.54 -4.54 -4.54 ...
$ slope : num 3.1 2.95 2.95 2.95 2.95 ...
$ female_mass : num 349 187 173 198 176 ...
$ hatchling_mass : num 100.9 69.2 73.7 81 84 ...
$ femass.log : num 1.7 0.8 0.7 0.9 0.7 0.9 0.5 1.6 1.4 1.3 ...
$ hatcmass.log : num 0.07 -0.42 -0.34 -0.23 -0.18 0.04 -0.45 0.74 0.24 -0.24 ...
$ rawFem : num 50.06 6.38 5.09 7.52 5.32 ...
$ rawHatc : num 1.17 0.38 0.45 0.59 0.65 1.1 0.35 5.44 1.74 0.57 ...
$ largest.clutch : int 17 8 8 8 7 8 6 1 2 18 ...
$ smallest.clutch : int 1 3 1 1 2 1 1 1 1 8 ...
$ largest.mean.clutch.size : num 11.4 5 6 4.4 4.8 3 3.3 1 2 14 ...
$ smallest.mean.clutch.size : num 11.4 4.8 2.6 3 2.5 2.6 2.5 1 2 11.5 ...
$ Longitude.centroid..from.Roll.et.al..2017.: num 30.3 34.8 19.8 -3.2 25.1 ...
$ Latitude.centroid..from.Roll.et.al..2017. : num -9.53 31.17 23.35 37.24 31.05 ...
$ Activity.time : chr "Diurnal" "Diurnal" "Diurnal" "Diurnal" ...
$ foraging.mode : chr "sit and wait" "widely foraging" "widely foraging" "mixed" ...
$ diet : chr "Carnivorous" "Carnivorous" "Carnivorous" "Carnivorous" ...
$ Binomial : chr "Acanthocercus_atricollis" "Acanthodactylus_beershebensis" "Acanthodactylus_boskianus" "Acanthodactylus_erythrurus" ...
$ M_females : num 22.2 15.3 15.4 15.3 15.4 ...
$ M_hatchlings : num 0.988 0.759 0.764 0.778 0.78 ...
$ mclutch : num 11.4 4.9 4.3 3.7 3.65 ...
$ rep_out : num 11.27 3.72 3.28 2.88 2.85 ...tapply(cdat$M_females, cdat$Family, length) Agamidae Anguidae Anniellidae Carphodactylidae
39 8 2 4
Chamaeleonidae Cordylidae Corytophanidae Crotaphytidae
5 6 3 5
Dactyloidae Diplodactylidae Eublepharidae Gekkonidae
40 14 6 35
Gerrhosauridae Gymnophthalmidae Helodermatidae Hoplocercidae
1 17 2 1
Lacertidae Leiocephalidae Liolaemidae Opluridae
49 3 11 1
Phrynosomatidae Phyllodactylidae Polychrotidae Pygopodidae
39 15 2 4
Scincidae Shinisauridae Sphaerodactylidae Teiidae
94 1 13 20
Tropiduridae Varanidae Xantusiidae Xenosauridae
12 28 3 2par(mar = c(1, 1, 1, 1))
plot(NA, NA, xlim = c(0, 10), ylim = c(-1, 10), ann = FALSE, axes = FALSE)
text(5, 10, "Foraging mode")
segments(3.5, 9.6, 6.5, 9.6, lwd = 2)
segments(3.5, 10.4, 6.5, 10.4, lwd = 2)
segments(3.5, 10.4, 3.5, 9.6, lwd = 2)
segments(6.5, 10.4, 6.5, 9.6, lwd = 2)
text(5, 7, "Distance traveled")
segments(3.5, 6.5, 6.5, 6.5, lwd = 2)
segments(3.5, 7.5, 6.5, 7.5, lwd = 2)
segments(3.5, 6.5, 3.5, 7.5, lwd = 2)
segments(6.5, 7.5, 6.5, 6.5, lwd = 2)
text(4.97, 3.7, "Energy gain")
segments(3.7, 4.2, 6.2, 4.2, lwd = 2)
segments(3.7, 3.2, 6.2, 3.2, lwd = 2)
segments(3.7, 4.2, 3.7, 3.2, lwd = 2)
segments(6.2, 4.2, 6.2, 3.2, lwd = 2)
text(8.6, 7, "Predation risk")
segments(7.5, 6.5, 9.8, 6.5, lwd = 2)
segments(7.5, 7.5, 9.8, 7.5, lwd = 2)
segments(7.5, 6.5, 7.5, 7.5, lwd = 2)
segments(9.8, 7.5, 9.8, 6.5, lwd = 2)
text(8.64, 3.7, "Reproductive output")
segments(7.2, 4.2, 10.1, 4.2, lwd = 2)
segments(7.2, 3.2, 10.1, 3.2, lwd = 2)
segments(7.2, 4.2, 7.2, 3.2, lwd = 2)
segments(10.1, 3.2, 10.1, 4.2, lwd = 2)
text(1.57, 3.7, "Energy expenditure")
segments(0, 4.2, 3.1, 4.2, lwd = 2)
segments(0, 3.2, 3.1, 3.2, lwd = 2)
segments(0, 4.2, 0, 3.2, lwd = 2)
segments(3.1, 4.2, 3.1, 3.2, lwd = 2)
text(4.9, 0, "Net energy")
segments(3.7, 0.5, 6.2, 0.5, lwd = 2)
segments(3.7, -0.5, 6.2, -0.5, lwd = 2)
segments(3.7, 0.5, 3.7, -0.5, lwd = 2)
segments(6.2, 0.5, 6.2, -0.5, lwd = 2)
Arrows(5, 9.5, 5, 7.8, lwd = 2, arr.type = "triangle")
Arrows(5, 6.3, 5, 4.8, lwd = 2, arr.type = "triangle")
Arrows(5, 4.8, 5, 6.1, lwd = 2, arr.type = "triangle")
Arrows(3.8, 6.3, 2, 4.8, lwd = 2, arr.type = "triangle")
Arrows(8.3, 4.76, 6.3, 6.2, lwd = 2, arr.type = "triangle")
Arrows(8.7, 6.1, 8.7, 4.8, lwd = 2, arr.type = "triangle")
Arrows(8.7, 4.8, 8.7, 6.1, lwd = 2, arr.type = "triangle")
Arrows(6.6, 7, 7, 7, lwd = 2, arr.type = "triangle")
Arrows(7.1, 3.7, 6.5, 3.7, lwd = 2, arr.type = "triangle")
Arrows(5, 3, 5, 1, lwd = 2, arr.type = "triangle")
Arrows(1.8, 3, 3.8, 1, lwd = 2, arr.type = "triangle")
Arrows(6.1, 1, 8.3, 3, lwd = 2, arr.type = "triangle")
Figure 1. A conceptual model depicting putative relationships among foraging behavior, energetics, predation risk, and reproductive output. The predicted relationships were derived from theoretical models of life-history evolution (see text for details).
layout(matrix(c(0, 1, 0,
0, 1, 0,
0, 2, 0,
0, 2, 0,
0, 3, 0,
0, 3, 0), nrow = 6, ncol = 3, byrow = TRUE))
## A
bertalanffy <- function(la, k, t){
la <- as.numeric(la)
k <- as.numeric(k)
t <- as.numeric(t)
la*(1 - 4*exp(-k*t))
}
fem1 <- bertalanffy(la = 4, k = 0.3, t = seq(1, 35, 1))
fem2 <- bertalanffy(la = 4.9, k = 0.25, t = seq(1, 35, 1))
par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0), cex.lab = 1.2)
plot(fem1 - 0.7 ~ seq(1, 35, 1), xlim = c(-2, 35), ylim = c(-6, 8), type = "l", lwd = 2, xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = "time (t)")
lines(seq(1, 35, 1) - 0.48, fem2, lwd = 2, lty = 4)
legend("topleft", legend = c("large female", "small female"), lty = c(4, 1), lwd = c(2, 2), bty = "n", cex = 0.7)
mtext("(a)", at = -6, line = 1, cex = 0.8, font = 2)
## B
fem3 <- bertalanffy(la = 5, k = 0.082, t = seq(1, 35, 1))
fem4 <- bertalanffy(la = 10, k = 0.21, t = seq(1, 35, 1))
par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0))
plot(fem3 - 0.7 ~ seq(1, 35, 1), xlim = c(4, 35), ylim = c(-6, 5), type = "l", lwd = 2, col = "skyblue", xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = " ")
lines(seq(1, 35, 1), fem4 - 8, lwd = 2, lty = 2, col = "purple")
segments(2.6, -6.5, 25, 6.2, lwd = 2)
segments(2.6, -6.5, 33.2, 6, lwd = 2)
segments(13.2, -0.5, 2.6, -0.5, lwd = 2, lty = 2, col = "purple")
segments(15.9, -1.1, 2.6, -1.1, lwd = 2, col = "skyblue")
segments(13.2, -0.5, 13.2, -7, lwd = 2, lty = 2, col = "purple")
segments(15.9, -1.1, 15.9, -7, lwd = 2, col = "skyblue")
legend("topleft", legend = c("sit and wait", "widely foraging"), lty = c(1, 2), lwd = c(2, 2), col = c("skyblue", "purple"), bty = "n", cex = 0.7)
mtext("(b)", at = 0.6, line = 1, cex = 0.8, font = 2)
mtext(expression("t"[wf]), side = 1, at = 13.2, line = 0.5, cex = 0.8, font = 2)
mtext(expression("t"[sw]), side = 1, at = 17.2, line = 0.5, cex = 0.8, font = 2)
mtext("time (t)", side = 1, at = 25, line = 0.9, cex = 0.8)
## C
par(mar = c(3, 3, 2, 1), mgp = c(0.8, 0, 0))
curve(2 - log(x)^2, ylim = c(-0.5, 3), xlim = c(0.1, 1), lwd = 2, lty = 2, col = "purple", xaxt = "n", yaxt = "n", ylab = "Cumulative energy gain (m)", xlab = "time (t)")
curve(log(x) + 1, add = TRUE, lwd = 2, col = "skyblue")
legend("topleft", legend = c("sit and wait", "widely foraging"), lty = c(1, 2), lwd = c(2, 2), col = c("skyblue", "purple"), bty = "n", cex = 0.7)
mtext("(c)", at = 0.001, line = 1, cex = 0.8, font = 2)
Figure 4. A theoretical model relating the cumulative energy gain for reproduction, m, as a function of time spent foraging, t, and maternal size. (a) Larger females reach their maximum capacity at a higher value of m than small females. (b) Widely-foraging females that are smaller but more efficient foragers may produce a greater reproductive output than larger sit-and-wait females. (c) Widely-foraging females may also produce a greater reproductive output than sit-and-wait females if they are both larger and more efficient foragers. twf and tsw in (b) represent the optimal foraging time of widely-foraging females and sit-and-wait females, respectively.
Fitting models (PGLS)
## Lacertilia
cdat$foraging.mode <- factor(cdat$foraging.mode, levels = c("widely foraging", "sit and wait", "mixed")) ## I am changing the order of the levels here so that the contrast of the coefficients can be made against the level "widely foraging".
rawMass <- cdat
rawMass <- rawMass[rawMass$rawFem < 5000, ] # Excluding a mixed-foraging species (too heavy)
check_rawMass <- name.check(tree, rawMass)
rm_phy_rawMass <- check_rawMass$tree_not_data
rm_dat_rawMass <- check_rawMass$data_not_tree
ctree_rawMass <- drop.tip(tree, rm_phy_rawMass)
name.check(ctree_rawMass, rawMass) [1] "OK"## Reproductive output
rawMass$rep_out_mass <- (rawMass$mclutch * rawMass$rawHatc)
## Number of individuals in each category (foraging mode)
tapply(rawMass$rawFem, rawMass$foraging.mode, length) widely foraging sit and wait mixed
202 222 60## Mean mass per category
tapply(rawMass$rawFem, rawMass$foraging.mode, mean) widely foraging sit and wait mixed
153.48926 32.35770 21.48317## Mean reproductive output per category
tapply(rawMass$rep_out_mass, rawMass$foraging.mode, mean) widely foraging sit and wait mixed
28.828361 5.540074 6.718067mod_mass <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corBrownian(phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass) Generalized least squares fit by maximum likelihood
Model: rep_out_mass ~ rawFem * foraging.mode
Data: rawMass
AIC BIC logLik
4938.271 4967.545 -2462.135
Correlation Structure: corBrownian
Formula: ~Binomial
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 6.963796 27.550624 0.252764 0.8006
rawFem 0.139790 0.005549 25.189779 0.0000
foraging.modesit and wait 0.809812 5.670050 0.142823 0.8865
foraging.modemixed -1.512868 6.179807 -0.244808 0.8067
rawFem:foraging.modesit and wait -0.081644 0.010811 -7.552212 0.0000
rawFem:foraging.modemixed 0.026696 0.089214 0.299233 0.7649
Correlation:
(Intr) rawFem frg.aw frgng. rF:.aw
rawFem -0.027
foraging.modesit and wait -0.105 0.099
foraging.modemixed -0.063 0.042 0.426
rawFem:foraging.modesit and wait 0.022 -0.308 -0.211 -0.088
rawFem:foraging.modemixed -0.004 0.001 0.106 -0.485 0.030
Standardized residuals:
Min Q1 Med Q3 Max
-2.21896992 -0.08372645 -0.07286486 -0.05356706 4.38147263
Residual standard error: 85.60208
Degrees of freedom: 484 total; 478 residualmod_mass1 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corBrownian(phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass1) Generalized least squares fit by maximum likelihood
Model: rep_out_mass ~ rawFem + foraging.mode
Data: rawMass
AIC BIC logLik
4989.078 5009.988 -2489.539
Correlation Structure: corBrownian
Formula: ~Binomial
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 11.599882 29.087362 0.398795 0.6902
rawFem 0.126833 0.005574 22.752688 0.0000
foraging.modesit and wait -8.569257 5.814017 -1.473896 0.1412
foraging.modemixed -4.064938 5.688570 -0.714580 0.4752
Correlation:
(Intr) rawFem frg.aw
rawFem -0.021
foraging.modesit and wait -0.102 0.036
foraging.modemixed -0.073 0.024 0.545
Standardized residuals:
Min Q1 Med Q3 Max
-2.67186452 -0.11700019 -0.04441283 -0.02384972 4.47349177
Residual standard error: 90.58862
Degrees of freedom: 484 total; 480 residualmod_mass2 <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corPagel(value = 0, phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
summary(mod_mass2) Generalized least squares fit by maximum likelihood
Model: rep_out_mass ~ rawFem * foraging.mode
Data: rawMass
AIC BIC logLik
4761.285 4794.742 -2372.642
Correlation Structure: corPagel
Formula: ~Binomial
Parameter estimate(s):
lambda
-0.08292133
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.221955 0.383885 0.57818 0.5634
rawFem 0.185596 0.002613 71.03444 0.0000
foraging.modesit and wait 1.803715 0.185985 9.69817 0.0000
foraging.modemixed -5.678696 3.737821 -1.51925 0.1294
rawFem:foraging.modesit and wait -0.122752 0.011515 -10.65994 0.0000
rawFem:foraging.modemixed 0.106016 0.112114 0.94561 0.3448
Correlation:
(Intr) rawFem frg.aw frgng. rF:.aw
rawFem -0.279
foraging.modesit and wait 0.168 0.096
foraging.modemixed -0.907 0.279 -0.482
rawFem:foraging.modesit and wait -0.006 -0.473 -0.560 0.021
rawFem:foraging.modemixed 0.042 -0.216 0.618 -0.425 0.022
Standardized residuals:
Min Q1 Med Q3 Max
-8.754409279 -0.041686637 -0.001960039 0.068491652 9.544845883
Residual standard error: 32.6597
Degrees of freedom: 484 total; 478 residualAnova(mod_mass2, type = "III") Analysis of Deviance Table (Type III tests)
Response: rep_out_mass
Df Chisq Pr(>Chisq)
(Intercept) 1 0.3343 0.5631
rawFem 1 5045.8917 <2e-16 ***
foraging.mode 2 107.0707 <2e-16 ***
rawFem:foraging.mode 2 115.0202 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_mass3 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corPagel(value = 0, phy = ctree_rawMass, form = ~Binomial), data = rawMass)
summary(mod_mass3) Generalized least squares fit by REML
Model: rep_out_mass ~ rawFem + foraging.mode
Data: rawMass
AIC BIC logLik
4871.726 4896.769 -2429.863
Correlation Structure: corPagel
Formula: ~Binomial
Parameter estimate(s):
lambda
-0.001644695
Coefficients:
Value Std.Error t-value p-value
(Intercept) 4.896045 2.608520 1.87694 0.0611
rawFem 0.155177 0.004527 34.27451 0.0000
foraging.modesit and wait -4.181294 3.548787 -1.17823 0.2393
foraging.modemixed -1.508168 5.441039 -0.27718 0.7818
Correlation:
(Intr) rawFem frg.aw
rawFem -0.263
foraging.modesit and wait -0.740 0.151
foraging.modemixed -0.492 0.110 0.348
Standardized residuals:
Min Q1 Med Q3 Max
-8.892572381 -0.119008564 -0.045067474 -0.002331544 9.525337598
Residual standard error: 36.92721
Degrees of freedom: 484 total; 480 residual## mod_mass2 is the most likely model even when the models below are included in the analysis
#mod_mass4 <- gls(rep_out_mass ~ rawFem*foraging.mode, correlation = corPagel(value = 1, phy = ctree_rawMass, form = ~Binomial), data = rawMass, method = "ML")
#summary(mod_mass4)
#mod_mass5 <- gls(rep_out_mass ~ rawFem + foraging.mode, correlation = corPagel(value = 1, phy = ctree_rawMass, form = ~Binomial), data = rawMass)
#summary(mod_mass5)
## Model selection table
output_rawMass <- model.sel(mod_mass, mod_mass1, mod_mass2, mod_mass3)
sel.table_mass <- round(as.data.frame(output_rawMass)[8:12], 3)
names(sel.table_mass)[1] <- "K"
sel.table_mass$Model <- rownames(sel.table_mass)
sel.table_mass$Model <- as.character(c(formula(mod_mass2), formula(mod_mass3), formula(mod_mass), formula(mod_mass1)))
sel.table_mass <- sel.table_mass[ , c(6, 1, 2, 3, 4, 5)]
sel.table_mass Model K logLik AICc delta
mod_mass2 rep_out_mass ~ rawFem * foraging.mode 8 -2372.642 4761.588 0.000
mod_mass3 rep_out_mass ~ rawFem + foraging.mode 6 -2429.863 4871.902 110.314
mod_mass rep_out_mass ~ rawFem * foraging.mode 7 -2462.135 4938.506 176.918
mod_mass1 rep_out_mass ~ rawFem + foraging.mode 5 -2489.539 4989.203 227.615
weight
mod_mass2 1
mod_mass3 0
mod_mass 0
mod_mass1 0Model checking (raw mass)
plot(mod_mass2, resid(., type = "n") ~ fitted(.), main = "Residuals v Fitted Values", abline = c(0, 0))
res <- resid(mod_mass2, type="n")
qqnorm(res, las = 1)
qqline(res)
A visual check of the most likely model indicates that the homogeneity and normality of residuals seem to be affected by influential cases (i.e., outliers). As a result, we standardized hatchling mass and maternal mass to alleviate this issue (see model checking below). In doing so, we recovered a more symetrical distibution of body mass among lizard species, increasing the robustness of our analysis. Influential cases should not be dropped from the model as there is no reasonable justification to do so. Instead, the patterns observed in our study could be resolved or even alter as more or better data become available.
## Standardized mass, M.
dat_full <- cdat
check_full <- name.check(tree, dat_full)
rm_phy_full <- check_full$tree_not_data
rm_dat_full <- check_full$data_not_tree
ctree_full <- drop.tip(tree, rm_phy_full)
name.check(ctree_full, dat_full) [1] "OK"## Number of individuals in each category (foraging mode)
tapply(dat_full$M_females, dat_full$foraging.mode, length) widely foraging sit and wait mixed
202 222 61## Mean scaled mass index per category
tapply(dat_full$M_females, dat_full$foraging.mode, mean) widely foraging sit and wait mixed
14.24854 17.88123 15.77951## Mean reproductive output per category
tapply(dat_full$rep_out, dat_full$foraging.mode, mean) widely foraging sit and wait mixed
3.220662 4.861416 5.788185full_mod <- gls(rep_out ~ M_females*foraging.mode, correlation = corBrownian(phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod) Generalized least squares fit by maximum likelihood
Model: rep_out ~ M_females * foraging.mode
Data: dat_full
AIC BIC logLik
2883.82 2913.109 -1434.91
Correlation Structure: corBrownian
Formula: ~Binomial
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 7.394834 3.950072 1.8720755 0.0618
M_females -0.222565 0.143833 -1.5473860 0.1224
foraging.modesit and wait 2.451992 1.976688 1.2404548 0.2154
foraging.modemixed -1.733346 2.026141 -0.8554912 0.3927
M_females:foraging.modesit and wait -0.197860 0.121228 -1.6321394 0.1033
M_females:foraging.modemixed 0.190587 0.126198 1.5102234 0.1316
Correlation:
(Intr) M_fmls frg.aw frgng. M_:.aw
M_females -0.557
foraging.modesit and wait -0.193 0.323
foraging.modemixed -0.141 0.266 0.442
M_females:foraging.modesit and wait 0.170 -0.339 -0.938 -0.500
M_females:foraging.modemixed 0.145 -0.309 -0.468 -0.947 0.587
Standardized residuals:
Min Q1 Med Q3 Max
-0.8510203 -0.3170531 -0.1719674 0.2257898 5.2981955
Residual standard error: 10.19489
Degrees of freedom: 485 total; 479 residualfull_mod1 <- gls(rep_out ~ M_females + foraging.mode, correlation = corBrownian(phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod1) Generalized least squares fit by maximum likelihood
Model: rep_out ~ M_females + foraging.mode
Data: dat_full
AIC BIC logLik
2891.793 2912.713 -1440.896
Correlation Structure: corBrownian
Formula: ~Binomial
Parameter estimate(s):
numeric(0)
Coefficients:
Value Std.Error t-value p-value
(Intercept) 7.826492 3.926231 1.993386 0.0468
M_females -0.242185 0.135280 -1.790247 0.0740
foraging.modesit and wait -1.189500 0.662400 -1.795743 0.0732
foraging.modemixed 1.604453 0.640735 2.504084 0.0126
Correlation:
(Intr) M_fmls frg.aw
M_females -0.537
foraging.modesit and wait -0.119 0.062
foraging.modemixed -0.036 -0.046 0.539
Standardized residuals:
Min Q1 Med Q3 Max
-0.8269391 -0.2711459 -0.1582826 0.2007470 5.4010162
Residual standard error: 10.3215
Degrees of freedom: 485 total; 481 residualfull_mod2 <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod2) Generalized least squares fit by maximum likelihood
Model: rep_out ~ M_females * foraging.mode
Data: dat_full
AIC BIC logLik
2770.838 2804.311 -1377.419
Correlation Structure: corPagel
Formula: ~Binomial
Parameter estimate(s):
lambda
0.8119803
Coefficients:
Value Std.Error t-value p-value
(Intercept) -0.194013 2.6555475 -0.0730595 0.9418
M_females 0.270669 0.1240541 2.1818666 0.0296
foraging.modesit and wait 3.101662 1.9832886 1.5638984 0.1185
foraging.modemixed -4.069565 2.2677650 -1.7945267 0.0734
M_females:foraging.modesit and wait -0.231656 0.1257082 -1.8428111 0.0660
M_females:foraging.modemixed 0.321048 0.1451586 2.2117076 0.0275
Correlation:
(Intr) M_fmls frg.aw frgng. M_:.aw
M_females -0.701
foraging.modesit and wait -0.328 0.452
foraging.modemixed -0.223 0.325 0.455
M_females:foraging.modesit and wait 0.311 -0.485 -0.949 -0.484
M_females:foraging.modemixed 0.221 -0.359 -0.467 -0.957 0.545
Standardized residuals:
Min Q1 Med Q3 Max
-1.2090603 -0.3604564 -0.2210199 0.1140136 7.4176637
Residual standard error: 6.419291
Degrees of freedom: 485 total; 479 residualfull_mod3 <- gls(rep_out ~ M_females + foraging.mode, correlation = corPagel(value = 0, phy = ctree_full, form = ~Binomial), data = dat_full, method = "ML")
summary(full_mod3) Generalized least squares fit by maximum likelihood
Model: rep_out ~ M_females + foraging.mode
Data: dat_full
AIC BIC logLik
2784.802 2809.907 -1386.401
Correlation Structure: corPagel
Formula: ~Binomial
Parameter estimate(s):
lambda
0.8188292
Coefficients:
Value Std.Error t-value p-value
(Intercept) 0.7964322 2.5818235 0.3084766 0.7579
M_females 0.2077960 0.1098458 1.8917060 0.0591
foraging.modesit and wait -0.8114631 0.6242388 -1.2999242 0.1942
foraging.modemixed 1.1099727 0.6600939 1.6815375 0.0933
Correlation:
(Intr) M_fmls frg.aw
M_females -0.657
foraging.modesit and wait -0.122 -0.006
foraging.modemixed -0.080 -0.010 0.513
Standardized residuals:
Min Q1 Med Q3 Max
-0.77190136 -0.34187016 -0.21073950 0.06296427 7.81257815
Residual standard error: 6.59119
Degrees of freedom: 485 total; 481 residual#full_mod4 <- gls(rep_out ~ M_females*foraging.mode, correlation = corPagel(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod4)
#full_mod5 <- gls(rep_out ~ M_females + foraging.mode, correlation = corPagel(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod5)
#full_mod6 <- gls(rep_out ~ M_females*foraging.mode, correlation = corMartins(value = 1, phy = ctree_full, form = ~Binomial, fixed = FALSE), data = dat_full, method = "ML")
#summary(full_mod6)Model checking (scaled mass index)
plot(full_mod2, resid(., type = "n") ~ fitted(.), main = "Residuals v Fitted Values", abline = c(0, 0))
res2 <- resid(full_mod2, type="n")
qqnorm(res2, las = 1)
qqline(res2)
Table S1. Model selection table of candidate models included in the analysis
output_full <- model.sel(full_mod, full_mod1, full_mod2, full_mod3)
sel.table_full <- round(as.data.frame(output_full)[7:11], 3)
names(sel.table_full)[1] <- "K"
sel.table_full$Model <- rownames(sel.table_full)
sel.table_full$Model <- as.character(c(formula(full_mod2), formula(full_mod3), formula(full_mod), formula(full_mod1)))
sel.table_full <- sel.table_full[ , c(6, 1, 2, 3, 4, 5)]
kable(sel.table_full, digits = 3, row.names = FALSE) %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")| Model | K | logLik | AICc | delta | weight |
|---|---|---|---|---|---|
| rep_out ~ M_females * foraging.mode | 8 | -1377.419 | 2771.141 | 0.000 | 0.999 |
| rep_out ~ M_females + foraging.mode | 6 | -1386.401 | 2784.977 | 13.837 | 0.001 |
| rep_out ~ M_females * foraging.mode | 7 | -1434.910 | 2884.055 | 112.915 | 0.000 |
| rep_out ~ M_females + foraging.mode | 5 | -1440.896 | 2891.918 | 120.777 | 0.000 |
Table 1. Analysis of Deviance Table (Type III tests) for the most likely model, based on the ranking of AICc for potential candidate models.
| Df | Chisq | Pr(>Chisq) | |
|---|---|---|---|
| (Intercept) | 1 | 0.005 | 0.942 |
| M_females | 1 | 4.761 | 0.029 |
| foraging.mode | 2 | 10.372 | 0.006 |
| M_females:foraging.mode | 2 | 18.105 | 0.000 |
layout(matrix(c(0, 0, 0, 0,
1, 1, 2, 2,
1, 1, 2, 2,
3, 3, 4, 4,
3, 3, 4, 4,
0, 0, 0, 0), nrow = 6, ncol = 4, byrow = TRUE))
par(mar = c(0, 5, 0, 0.5))
## Distribution of body size
rm1 <- density(rawMass$rawFem[rawMass$foraging.mode == "widely foraging"])
rm2 <- density(rawMass$rawFem[rawMass$foraging.mode == "mixed"])
rm3 <- density(rawMass$rawFem[rawMass$foraging.mode == "sit and wait"])
plot(1, type = "n", ann = FALSE, bty = "n", xaxt = "n", yaxt = "n")
par(mar = c(0, 4, 0, 0.5))
t1 <- density(dat_full$M_females[dat_full$foraging.mode == "widely foraging"])
t2 <- density(dat_full$M_females[dat_full$foraging.mode == "mixed"])
t3 <- density(dat_full$M_females[dat_full$foraging.mode == "sit and wait"])
plot(t1, ylim = c(0, max(c(t1$y, t1$y))), ann = FALSE, bty = "n", xaxt = "n", yaxt = "n")
polygon(t1, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(t2)
polygon(t2, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(t3)
polygon(t3, col = rgb(0, 0, 1, alpha = 0.2))
## Lacertilia
## Raw mass
## widely foraging
rw_fg <- rawMass[rawMass$foraging.mode == "widely foraging", ]
check_rm_wf <- name.check(tree, rw_fg)
rm_phy_rm_wf <- check_rm_wf$tree_not_data
rm_dat_rm_wf <- check_rm_wf$data_not_tree
ctree_rm_wf <- drop.tip(tree, rm_phy_rm_wf)
#name.check(ctree_rm_wf, rw_fg)
SSX.rm_wf <- sum(round((rw_fg$rawFem - mean(rw_fg$rawFem))^2), 2)
s2.rm_wf <- var(rw_fg$rep_out_mass)
n.rm_wf <- length(rw_fg$rep_out_mass)
x.rm_wf <- seq(0.50, 4072.0, 1)
m.x.rm_wf <- mean(round(rw_fg$rawFem,1))
se.rm_wf <- sqrt(s2.rm_wf*((1/n.rm_wf) + (((x.rm_wf - m.x.rm_wf)^2)/SSX.rm_wf)))
is.rm_wf <- qt(0.975, df = n.rm_wf - 2)
ii.rm_wf <- qt(0.025, df = n.rm_wf - 2)
ic.s.rm_wf <- se.rm_wf*is.rm_wf
ic.i.rm_wf <- se.rm_wf*ii.rm_wf
upper.i.rm_wf <- (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf) + ic.s.rm_wf
lower.i.rm_wf <- (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf) + ic.i.rm_wf
# sit and wait
rm_sw <- rawMass[rawMass$foraging.mode == "sit and wait", ]
check_rm_sw <- name.check(tree, rm_sw)
rm_phy_rm_sw <- check_rm_sw$tree_not_data
rm_dat_rm_sw <- check_rm_sw$data_not_tree
ctree_rm_sw <- drop.tip(tree, rm_phy_rm_sw)
#name.check(ctree_rm_sw, rm_sw)
SSX.rm_sw <- sum(round((rm_sw$rawFem - mean(rm_sw$rawFem))^2), 2)
s2.rm_sw <- var(rm_sw$rep_out_mass)
n.rm_sw <- length(rm_sw$rep_out_mass)
x.rm_sw <- seq(0.28, 3127.7, 1)
m.x.rm_sw <- mean(round(rm_sw$rawFem, 1))
se.rm_sw <- sqrt(s2.rm_sw*((1/n.rm_sw) + (((x.rm_sw - m.x.rm_sw)^2)/SSX.rm_sw)))
is.rm_sw <- qt(0.975, df = n.rm_sw - 2)
ii.rm_sw <- qt(0.025, df = n.rm_sw - 2)
ic.s.rm_sw <- se.rm_sw*is.rm_sw
ic.i.rm_sw <- se.rm_sw*ii.rm_sw
upper.i.rm_sw <- ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw) + ic.s.rm_sw
lower.i.rm_sw <- ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw) + ic.i.rm_sw
## mixed
rm_mx <- rawMass[rawMass$foraging.mode == "mixed", ]
check_rm_mx <- name.check(tree, rm_mx)
rm_phy_rm_mx <- check_rm_mx$tree_not_data
rm_dat_rm_mx <- check_rm_mx$data_not_tree
ctree_rm_mx <- drop.tip(tree, rm_phy_rm_mx)
#name.check(ctree_rm_mx, rm_mx)
SSX.rm_mx <- sum(round((rm_mx$rawFem - mean(rm_mx$rawFem))^2), 2)
s2.rm_mx <- var(rm_mx$rep_out_mass)
n.rm_mx <- length(rm_mx$rep_out_mass)
x.rm_mx <- seq(0.44, 222.1, 1)
m.x.rm_mx <- mean(round(rm_mx$rawFem, 1))
se.rm_mx <- sqrt(s2.rm_mx*((1/n.rm_mx) + (((x.rm_mx - m.x.rm_mx)^2)/SSX.rm_mx)))
is.rm_mx <- qt(0.975, df = n.rm_mx - 2)
ii.rm_mx <- qt(0.025, df = n.rm_mx - 2)
ic.s.rm_mx <- se.rm_mx*is.rm_mx
ic.i.rm_mx <- se.rm_mx*ii.rm_mx
upper.i.rm_mx.f <- ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx) + ic.s.rm_mx
lower.i.rm_mx.f <- ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx) + ic.i.rm_mx
par(mar = c(0, 5, 0, 0.5), mgp = c(3.2, 1, 0), las = 1)
plot(rawMass$rep_out_mass ~ rawMass$rawFem, type = "p", xlab = "mass (g)", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.5, pch = 21, bg = c("purple", "skyblue", "red")[as.numeric(rawMass$foraging.mode)], col = c("purple", "skyblue", "red")[as.numeric(rawMass$foraging.mode)])
polygon(c(rev(x.rm_mx), x.rm_mx), c(rev(lower.i.rm_mx.f), upper.i.rm_mx.f), border = FALSE, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(x = x.rm_mx, y = ((coef(mod_mass2)[1] + coef(mod_mass2)[4]) + (coef(mod_mass2)[2] + coef(mod_mass2)[6])*x.rm_mx), col = "red", lwd = 2, lty = 15)
polygon(c(rev(x.rm_sw), x.rm_sw), c(rev(lower.i.rm_sw), upper.i.rm_sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
lines(x = x.rm_sw, y = ((coef(mod_mass2)[1] + coef(mod_mass2)[3]) + (coef(mod_mass2)[2] + coef(mod_mass2)[5])*x.rm_sw), col = "skyblue", lwd = 2, lty = 2)
polygon(c(rev(x.rm_wf), x.rm_wf), c(rev(lower.i.rm_wf), upper.i.rm_wf), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(x = x.rm_wf, y = (coef(mod_mass2)[1] + coef(mod_mass2)[2]*x.rm_wf), col = "purple", lwd = 2, lty = 2)
mtext("(a)", side = 3, line = 1, at = -600, font = 2)
mtext("mass (g)", side = 1, line = 4, cex = 1)
legend("topleft", legend = c("mixed", "sit and wait", "widely foraging"), lty = c(15, 1, 2), bty = "n", pch = c(16, 16, 16), bg = c("black", "black", "black"), col = c("red", "skyblue", "purple"), cex = 1)
## Scaled mass index
## widely foraging
wf <- dat_full[dat_full$foraging.mode == "widely foraging", ]
check_wf <- name.check(tree, wf)
rm_phy_wf <- check_wf$tree_not_data
rm_dat_wf <- check_wf$data_not_tree
ctree_wf <- drop.tip(tree, rm_phy_wf)
#name.check(ctree_wf, wf)
SSX.wf <- sum(round((wf$M_females - mean(wf$M_females))^2), 2)
s2.wf <- var(wf$rep_out)
n.wf <- length(wf$rep_out)
x.wf <- seq(1.2, 31.8, 1)
m.x.wf <- mean(round(wf$M_females,1))
se.wf <- sqrt(s2.wf*((1/n.wf) + (((x.wf - m.x.wf)^2)/SSX.wf)))
is.wf <- qt(0.975, df = n.wf - 2)
ii.wf <- qt(0.025, df = n.wf - 2)
ic.s.wf <- se.wf*is.wf
ic.i.wf <- se.wf*ii.wf
upper.i <- (coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf) + ic.s.wf
lower.i <- (coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf) + ic.i.wf
## sit and wait
sw <- dat_full[dat_full$foraging.mode == "sit and wait", ]
check_sw <- name.check(tree, sw)
rm_phy_sw <- check_sw$tree_not_data
rm_dat_sw <- check_sw$data_not_tree
ctree_sw <- drop.tip(tree, rm_phy_sw)
#name.check(ctree_sw, sw)
SSX.sw <- sum(round((sw$M_females - mean(sw$M_females))^2), 2)
s2.sw <- var(sw$rep_out)
n.sw <- length(sw$rep_out)
x.sw <- seq(1.3, 30.8, 1)
m.x.sw <- mean(round(sw$M_females, 1))
se.sw <- sqrt(s2.sw*((1/n.sw) + (((x.sw - m.x.sw)^2)/SSX.sw)))
is.sw <- qt(0.975, df = n.sw - 2)
ii.sw <- qt(0.025, df = n.sw - 2)
ic.s.sw <- se.sw*is.sw
ic.i.sw <- se.sw*ii.sw
upper.i.sw <- ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw) + ic.s.sw
lower.i.sw <- ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw) + ic.i.sw
## mixed
mx <- dat_full[dat_full$foraging.mode == "mixed", ]
check_mx <- name.check(tree, mx)
rm_phy_mx <- check_mx$tree_not_data
rm_dat_mx <- check_mx$data_not_tree
ctree_mx <- drop.tip(tree, rm_phy_mx)
#name.check(ctree_mx, mx)
SSX.mx <- sum(round((mx$M_females - mean(mx$M_females))^2), 2)
s2.mx <- var(mx$rep_out)
n.mx <- length(mx$rep_out)
x.mx <- seq(1.0, 27.7, 1)
m.x.mx <- mean(round(mx$M_females, 1))
se.mx <- sqrt(s2.mx*((1/n.mx) + (((x.mx - m.x.mx)^2)/SSX.mx)))
is.mx <- qt(0.975, df = n.mx - 2)
ii.mx <- qt(0.025, df = n.mx - 2)
ic.s.mx <- se.mx*is.mx
ic.i.mx <- se.mx*ii.mx
upper.i.mx.f <- ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx) + ic.s.mx
lower.i.mx.f <- ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx) + ic.i.mx
par(mar = c(0, 4, 0, 0.5), mgp = c(2.8, 1, 0), las = 1)
plot(dat_full$rep_out ~ dat_full$M_females, type = "p", xlab = " ", ylab = "Hatchling mass x clutch size (g)", cex.lab = 1.5, pch = 21, bg = c("purple", "skyblue", "red")[as.numeric(dat_full$foraging.mode)], col = c("purple", "skyblue", "red")[as.numeric(dat_full$foraging.mode)])
polygon(c(rev(x.mx), x.mx), c(rev(lower.i.mx.f), upper.i.mx.f), border = FALSE, col = rgb(0.8, 0, 0, alpha = 0.3))
lines(x = x.mx, y = ((coef(full_mod2)[1] + coef(full_mod2)[4]) + (coef(full_mod2)[2] + coef(full_mod2)[6])*x.mx), col = "red", lwd = 2, lty = 15)
polygon(c(rev(x.sw), x.sw), c(rev(lower.i.sw), upper.i.sw), border = FALSE, col = rgb(0, 0, 1, alpha = 0.2))
lines(x = x.sw, y = ((coef(full_mod2)[1] + coef(full_mod2)[3]) + (coef(full_mod2)[2] + coef(full_mod2)[5])*x.sw), col = "skyblue", lwd = 2, lty = 2)
polygon(c(rev(x.wf), x.wf), c(rev(lower.i), upper.i), border = FALSE, col = rgb(0.8, 0, 0.7, alpha = 0.3))
lines(x = x.wf, y = coef(full_mod2)[1] + coef(full_mod2)[2]*x.wf, col = "purple", lwd = 2, lty = 2)
mtext("(b)", side = 3, line = 1, at = -2, font = 2)
mtext(expression(hat(M) (g)), side = 1, line = 4, cex = 1)
legend("topleft", legend = c("mixed", "sit-and-wait foraging", "widely foraging"), lty = c(15, 1, 2), bty = "n", pch = c(16, 16, 16), bg = c("black", "black", "black"), col = c("red", "skyblue", "purple"), cex = 1)
Figure 3. Effects of maternal mass and foraging mode on the evolution of reproductive output of lizards, as determined by phylogenetic generalized least squares analysis. (a) Difference in reproductive output among foraging modes assuming body mass as a proxy for energy. (b) Difference in reproductive output among foraging modes assuming the scaled mass index as a proxy for energy.
Ancenstral state reconstruction
dev.off() null device
1set.seed(1094)
## Ancestral character state
anc_st <- cdat[ , c(25, 23)]
anc_st[486, ] <- c("Sphenodon_punctatus", "sit and wait")
rownames(anc_st)[486] <- "Sphenodon_punctatus"
check_st <- name.check(tree, anc_st)
rm_phy_st <- check_st$tree_not_data
rm_dat_st <- check_st$data_not_tree
ctree_st <- drop.tip(tree, rm_phy_st)
anc_st <- subset(anc_st, subset = rownames(anc_st) %in% ctree_st$tip, select = names(anc_st)[1:2])
#name.check(ctree, anc_st)
forg <- as.matrix(anc_st)[,2]
head(forg) Acanthocercus_atricollis Acanthodactylus_beershebensis
"sit and wait" "widely foraging"
Acanthodactylus_boskianus Acanthodactylus_erythrurus
"widely foraging" "mixed"
Acanthodactylus_pardalis Acanthodactylus_schreiberi
"sit and wait" "widely foraging"## Reproductive output
anc_st_rp <- cdat[ , c(29, 23)]
anc_st_rp[486, 1] <- 0
anc_st_rp[486, 2] <- "sit and wait"
rownames(anc_st_rp)[486] <- "Sphenodon_punctatus"
check_rp <- name.check(tree, anc_st_rp)
rm_phy_rp <- check_rp$tree_not_data
rm_dat_rp <- check_rp$data_not_tree
ctree_rp <- drop.tip(tree, rm_phy_rp)
anc_st_rp <- subset(anc_st_rp, subset = rownames(anc_st_rp) %in% ctree_rp$tip, select = names(anc_st_rp)[1:2])
#name.check(ctree, anc_st_rp)
anc_st_rp$rep_out <- as.numeric(anc_st_rp$rep_out)
anc_st_rp$foraging.mode <- as.factor(anc_st_rp$foraging.mode)
anc_st_rp_bars <- anc_st_rp[1]
anc_st_rp_bars <- as.matrix(anc_st_rp_bars)[ , 1]
head(anc_st_rp_bars) Acanthocercus_atricollis Acanthodactylus_beershebensis
11.268193 3.720936
Acanthodactylus_boskianus Acanthodactylus_erythrurus
3.284710 2.879015
Acanthodactylus_pardalis Acanthodactylus_schreiberi
2.845583 2.252075## Reproductive output for widely foraging and sit-and-wait species
anc_st_rp2 <- anc_st_rp[anc_st_rp$foraging.mode != "mixed", ]
check_rp2 <- name.check(tree, anc_st_rp2)
rm_phy_rp2 <- check_rp2$tree_not_data
rm_dat_rp2 <- check_rp2$data_not_tree
ctree_rp2 <- drop.tip(tree, rm_phy_rp2)
anc_st_rp2 <- subset(anc_st_rp2, subset = rownames(anc_st_rp2) %in% ctree_rp2$tip, select = names(anc_st_rp2)[1:2])
#name.check(ctree_rp2, anc_st_rp2)
anc_st_rp2$rep_out <- as.numeric(anc_st_rp2$rep_out)
anc_st_rp2$foraging.mode <- as.factor(anc_st_rp2$foraging.mode)
anc_st_rp2_bars <- anc_st_rp2[1]
anc_st_rp2_bars <- as.matrix(anc_st_rp2_bars)[ , 1]
head(anc_st_rp2_bars) Acanthocercus_atricollis Acanthodactylus_beershebensis
11.2681934 3.7209360
Acanthodactylus_boskianus Acanthodactylus_pardalis
3.2847103 2.8455830
Acanthodactylus_schreiberi Acontias_gariepensis
2.2520747 0.7391738## Use of function ace (for comparison)
fitER <- ace(forg, ctree_st, model = "ER", type = "discrete")
fitARD <- ace(forg, ctree_st, model = "ARD", type = "discrete")
fitSR <- ace(forg, ctree_st, model = "SYM", type = "discrete")
fitER
Ancestral Character Estimation
Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "ER")
Log-likelihood: -335.5832
Rate index matrix:
mixed sit and wait widely foraging
mixed . 1 1
sit and wait 1 . 1
widely foraging 1 1 .
Parameter estimates:
rate index estimate std-err
1 0.0032 3e-04
Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
mixed sit and wait widely foraging
0.1631016 0.4946075 0.3422909fitARD
Ancestral Character Estimation
Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "ARD")
Log-likelihood: -320.1167
Rate index matrix:
mixed sit and wait widely foraging
mixed . 3 5
sit and wait 1 . 6
widely foraging 2 4 .
Parameter estimates:
rate index estimate std-err
1 0.0062 0.0011
2 0.0045 0.0010
3 0.0133 0.0040
4 0.0022 0.0007
5 0.0205 0.0051
6 0.0000 0.0007
Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
mixed sit and wait widely foraging
0.3347980 0.3811181 0.2840839fitSR
Ancestral Character Estimation
Call: ace(x = forg, phy = ctree_st, type = "discrete", model = "SYM")
Log-likelihood: -333.6765
Rate index matrix:
mixed sit and wait widely foraging
mixed . 1 2
sit and wait 1 . 3
widely foraging 2 3 .
Parameter estimates:
rate index estimate std-err
1 0.0041 6e-04
2 0.0041 7e-04
3 0.0023 4e-04
Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
mixed sit and wait widely foraging
0.2145978 0.4773200 0.3080822head(round(fitER$lik.anc, 3)) mixed sit and wait widely foraging
487 0.163 0.495 0.342
488 0.018 0.582 0.399
489 0.005 0.954 0.041
490 0.010 0.969 0.021
491 0.017 0.953 0.030
492 0.069 0.921 0.011head(round(fitARD$lik.anc, 3)) mixed sit and wait widely foraging
487 0.335 0.381 0.284
488 0.283 0.362 0.356
489 0.074 0.895 0.031
490 0.093 0.893 0.013
491 0.132 0.853 0.015
492 0.138 0.855 0.007head(round(fitSR$lik.anc, 3)) mixed sit and wait widely foraging
487 0.215 0.477 0.308
488 0.060 0.566 0.374
489 0.017 0.953 0.030
490 0.024 0.962 0.014
491 0.036 0.943 0.021
492 0.089 0.904 0.007Table S2.
ic <- data.frame(Model = c("ER", "ARD", "SYM"), AIC = c(AIC(fitER), AIC(fitARD), AIC(fitSR)), stringsAsFactors = FALSE)
kable(ic, digits = 3, row.names = FALSE) %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")| Model | AIC |
|---|---|
| ER | 673.166 |
| ARD | 652.233 |
| SYM | 673.353 |
## Overlaying posterior probabilities on the tree. This is done for the best fit model.
colors <- c("red", "skyblue", "purple")
cols <- setNames(colors[1:length(unique(forg))], sort(unique(forg)))
plotTree(ctree_st, type = "fan", fsize = 0.2, ftype = "i", part = 0.98)
t <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(115, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t - scale, cex = 0.6)
text(mean(scale), -38, "time (mya)", cex = 0.7)
nodelabels(node = 1:ctree_st$Nnode + Ntip(ctree_st), pie = fitARD$lik.anc, piecol = cols, cex = 0.2)
tiplabels(pie = to.matrix(forg[ctree_st$tip.label], unique(sort(forg))), piecol = cols, cex = 0.1)
add.simmap.legend(colors = cols, prompt = FALSE, x= 0.9*par()$usr[1], y = -max(nodeHeights(ctree_st)), fsize = 0.8)
ASER <- make.simmap(ctree_st, forg, model = "ER", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
mixed sit and wait widely foraging
mixed -0.006716348 0.003358174 0.003358174
sit and wait 0.003358174 -0.006716348 0.003358174
widely foraging 0.003358174 0.003358174 -0.006716348
and (mean) root node prior probabilities
pi =
mixed sit and wait widely foraging
0.3333333 0.3333333 0.3333333ASARD <- make.simmap(ctree_st, forg, model = "ARD", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
mixed sit and wait widely foraging
mixed -4.2957510 3.8697964 0.4259546
sit and wait 0.6786649 -2.1882412 1.5095763
widely foraging 0.5975036 0.9830754 -1.5805790
and (mean) root node prior probabilities
pi =
mixed sit and wait widely foraging
0.1291999 0.4275951 0.4432051ASSR <- make.simmap(ctree_st, forg, model = "SYM", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
mixed sit and wait widely foraging
mixed -1.56999075 0.01032159 1.5596692
sit and wait 0.01032159 -0.02980089 0.0194793
widely foraging 1.55966916 0.01947930 -1.5791485
and (mean) root node prior probabilities
pi =
mixed sit and wait widely foraging
0.3333333 0.3333333 0.3333333Table S3.
as <- data.frame(Model = c("modER", "modARD", "modSYM"), AIC = c(AIC(ASER), AIC(ASARD), AIC(ASSR)), stringsAsFactors = FALSE)
kable(as, digits = 3, row.names = FALSE) %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")| Model | AIC |
|---|---|
| modER | 685.665 |
| modARD | 969.306 |
| modSYM | 947.630 |
## Simulate 1000 stochastic character maps using empirical Bayes method
ASER <- make.simmap(ctree_st, forg, model = "ER", nsim = 1000, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 1e+05 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
mixed sit and wait widely foraging
mixed -0.006463354 0.003231677 0.003231677
sit and wait 0.003231677 -0.006463354 0.003231677
widely foraging 0.003231677 0.003231677 -0.006463354
and (mean) root node prior probabilities
pi =
mixed sit and wait widely foraging
0.3333333 0.3333333 0.3333333pd <- summary(ASER, plot = FALSE)
pd 1000 trees with a mapped discrete character with states:
mixed, sit and wait, widely foraging
trees have 126.77 changes between states on average
changes are of the following types:
mixed,sit and wait mixed,widely foraging sit and wait,mixed
x->y 8.05 8.06 34.261
sit and wait,widely foraging widely foraging,mixed
x->y 24.335 25.801
widely foraging,sit and wait
x->y 26.263
mean total time spent in each state is:
mixed sit and wait widely foraging total
raw 1.648428e+03 1.033519e+04 7900.6573103 19884.28
prop 8.290106e-02 5.197671e-01 0.3973319 1.00save <- countSimmap(ASER)
min(save$Tr[,1]) [1] 105write.table(pd$ace, "likelihood_nodes_ASER.csv", sep = ",")
## What happens if the tuatara is excluded?
forgt <- forg[-486]
ASERT <- make.simmap(ctree, forgt, model = "ER", nsim = 1000, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 1e+05 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
mixed sit and wait widely foraging
mixed -0.006425791 0.003212896 0.003212896
sit and wait 0.003212896 -0.006425791 0.003212896
widely foraging 0.003212896 0.003212896 -0.006425791
and (mean) root node prior probabilities
pi =
mixed sit and wait widely foraging
0.3333333 0.3333333 0.3333333pdt <- summary(ASERT, plot = FALSE)
pdt 1000 trees with a mapped discrete character with states:
mixed, sit and wait, widely foraging
trees have 124.251 changes between states on average
changes are of the following types:
mixed,sit and wait mixed,widely foraging sit and wait,mixed
x->y 7.342 7.678 33.895
sit and wait,widely foraging widely foraging,mixed
x->y 24.041 25.419
widely foraging,sit and wait
x->y 25.876
mean total time spent in each state is:
mixed sit and wait widely foraging total
raw 1.588090e+03 1.012479e+04 7818.285105 19531.16
prop 8.131054e-02 5.183915e-01 0.400298 1.00savet <- countSimmap(ASERT)
min(savet$Tr[,1]) [1] 102write.table(pdt$ace, "likelihood_nodes_ASERT.csv", sep = ",")
## Phylogenetic signal of foraging modes
anc_st_rp2$Binomial <- rownames(anc_st_rp2)
phylo.d(anc_st_rp2, ctree_rp2, Binomial, binvar = foraging.mode, permut = 5000)
Calculation of D statistic for the phylogenetic structure of a binary variable
Data : data
Binary variable : foraging.mode
Counts of states: widely foraging = 202
sit and wait = 223
Phylogeny : phy
Number of permutations : 5000
Estimated D : -0.144269
Probability of E(D) resulting from no (random) phylogenetic structure : 0
Probability of E(D) resulting from Brownian phylogenetic structure : 0.8004## Estimating ancestral character state only for widely foraging and sit-and-wait species
forg2 <- as.matrix(anc_st_rp2)[ , 2]
head(forg2) Acanthocercus_atricollis Acanthodactylus_beershebensis
"sit and wait" "widely foraging"
Acanthodactylus_boskianus Acanthodactylus_pardalis
"widely foraging" "sit and wait"
Acanthodactylus_schreiberi Acontias_gariepensis
"widely foraging" "widely foraging"ASER2 <- make.simmap(ctree_rp2, forg2, model = "ER", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
sit and wait widely foraging
sit and wait -0.002167574 0.002167574
widely foraging 0.002167574 -0.002167574
and (mean) root node prior probabilities
pi =
sit and wait widely foraging
0.5 0.5ASARD2 <- make.simmap(ctree_rp2, forg2, model = "ARD", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
sit and wait widely foraging
sit and wait -0.002013939 0.002013939
widely foraging 0.007239599 -0.007239599
and (mean) root node prior probabilities
pi =
sit and wait widely foraging
0.7823601 0.2176399ASSR2 <- make.simmap(ctree_rp2, forg2, model = "SYM", nsim = 1, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 100 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
sit and wait widely foraging
sit and wait -0.00306315 0.00306315
widely foraging 0.00306315 -0.00306315
and (mean) root node prior probabilities
pi =
sit and wait widely foraging
0.5 0.5Table S4.
as2 <- data.frame(Model = c("modER2", "modARD2", "modSR2"), AIC = c(AIC(ASER2), AIC(ASARD2), AIC(ASSR2)), stringsAsFactors = FALSE)
kable(as2, digits = 3, row.names = FALSE) %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")| Model | AIC |
|---|---|
| modER2 | 297.869 |
| modARD2 | 309.341 |
| modSR2 | 295.076 |
## simulate 1000 stochastic character maps using empirical Bayes method (excluding mixed foragers)
ASSYM2 <- make.simmap(ctree_rp2, forg2, model = "SYM", nsim = 1000, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 1e+05 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
sit and wait widely foraging
sit and wait -0.003025366 0.003025366
widely foraging 0.003025366 -0.003025366
and (mean) root node prior probabilities
pi =
sit and wait widely foraging
0.5 0.5pd2 <- summary(ASSYM2, plot = FALSE)
pd2 1000 trees with a mapped discrete character with states:
sit and wait, widely foraging
trees have 52.577 changes between states on average
changes are of the following types:
sit and wait,widely foraging widely foraging,sit and wait
x->y 26.137 26.44
mean total time spent in each state is:
sit and wait widely foraging total
raw 1.010405e+04 7712.4930593 17816.54
prop 5.671161e-01 0.4328839 1.00save2 <- countSimmap(ASSYM2)
min(save2$Tr[,1]) [1] 42write.table(pd2$ace, "likelihood_nodes_ASSYM2.csv", sep = ",")
## simulate 1000 stochastic character maps using empirical Bayes method (excluding mixed foragers and the tuatara)
forg2T <- cdat[ , c(25, 23)]
forg2T <- forg2T[forg2T$foraging.mode != "mixed", ]
check_st2 <- name.check(tree, forg2T)
rm_phy_st2 <- check_st2$tree_not_data
rm_dat_st2 <- check_st2$data_not_tree
ctree_st2 <- drop.tip(tree, rm_phy_st2)
forg2T <- subset(forg2T, subset = rownames(forg2T) %in% ctree_st2$tip, select = names(forg2T)[1:2])
#name.check(ctree_st2, forg2T)
forg2T <- as.matrix(forg2T)[,2]
head(forg2T) Acanthocercus_atricollis Acanthodactylus_beershebensis
"sit and wait" "widely foraging"
Acanthodactylus_boskianus Acanthodactylus_pardalis
"widely foraging" "sit and wait"
Acanthodactylus_schreiberi Acontias_gariepensis
"widely foraging" "widely foraging"ASSYMT <- make.simmap(ctree_st2, forg2T, model = "SYM", nsim = 1000, pi = "estimated", Q = "mcmc") Running MCMC burn-in. Please wait....
Running 1e+05 generations of MCMC, sampling every 100 generations.
Please wait....
make.simmap is simulating with a sample of Q from
the posterior distribution
Mean Q from the posterior is
Q =
sit and wait widely foraging
sit and wait -0.002965521 0.002965521
widely foraging 0.002965521 -0.002965521
and (mean) root node prior probabilities
pi =
sit and wait widely foraging
0.5 0.5pd2T <- summary(ASSYMT, plot = FALSE)
pd2T 1000 trees with a mapped discrete character with states:
sit and wait, widely foraging
trees have 51.071 changes between states on average
changes are of the following types:
sit and wait,widely foraging widely foraging,sit and wait
x->y 24.979 26.092
mean total time spent in each state is:
sit and wait widely foraging total
raw 9788.337854 7675.089625 17463.43
prop 0.560505 0.439495 1.00save2T <- countSimmap(ASSYMT)
min(save2T$Tr[,1]) [1] 42write.table(pd2T$ace, "likelihood_nodes_ASSYMT.csv", sep = ",")
## Plotting trees
colors <- c("red", "skyblue", "purple")
cols <- setNames(colors[1:length(unique(forg))], sort(unique(forg)))
colorbars <-setNames(c("red", "skyblue", "purple"), c("mixed", "sit and wait", "widely foraging"))
ss <-getStates(ASER, "tips")
barcols <- setNames(sapply(ss, function(x, y) y[which(names(y) == x)], y = colorbars), names(ss))
i <- 1:1000
plotTree.wBars(ASER[[sample(i, 1)]], log1p(anc_st_rp_bars), type = "fan", method = "plotSimmap", colors = cols, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols, border = FALSE)
t1 <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t1, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(119, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t1 - scale, cex = 0.9)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd$ace, piecol = cols, cex = 0.26)
add.simmap.legend(colors = cols, prompt = FALSE, x = -120, y = 140, fsize = 0.9)
Figure 2. Random sample of stochastic character maps depicting the evolution of foraging mode in 485 species of lizards. Bars at the tips of the phylogeny represent log-transformed values of reproductive output for all lizards, but not the outgroup,Sphenodon punctatus. Pie charts on internal nodes represent posterior probability estimates resulting from 1,000 simulations of the character histories. Major clades are enumerated as follows: 1) Gekkota, 2) Scincoidea, 3) Lacertoidea, 4) Anguimorpha, 5) Toxicofera, and 6) Iguania. Lizard photos by Mark O’Shea.
## Ancestral state reconstruction excluding the tuatara
colst <- setNames(colors[1:length(unique(forgt))], sort(unique(forgt)))
sst <-getStates(ASERT, "tips")
colorbars <-setNames(c("red", "skyblue", "purple"), c("mixed", "sit and wait", "widely foraging"))
barcolst <- setNames(sapply(sst, function(x, y) y[which(names(y) == x)], y = colorbars), names(sst))
iT <- 1:1000
plotTree.wBars(ASERT[[sample(iT, 1)]], log1p(anc_st_rp_bars[-486]), type = "fan", method = "plotSimmap", colors = colst, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcolst, border = FALSE)
t2 <- max(nodeHeights(ctree))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t2, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(55, 20, "1", cex = 1)
text(10, 20, "2", cex = 1)
text(-30, 8, "3", cex = 1)
text(-72, -38, "4", cex = 1)
text(-11, -24, "5", cex = 1)
text(-8, -40, "6", cex = 1)
text(scale, rep(-20, length(scale)), t2 - scale, cex = 0.7)
text(mean(scale), -35, "time (mya)", cex = 0.9)
nodelabels(pie = pdt$ace, piecol = colst, cex = 0.3)
add.simmap.legend(colors = colst, prompt = FALSE, x = -80, y = 90, fsize = 0.9)
## Internal nodes numbering
plotTree.wBars(ASER[[sample(i, 1)]], log1p(anc_st_rp_bars), type = "fan", method = "plotSimmap", colors = cols, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols, border = FALSE)
t1 <- max(nodeHeights(ctree_st))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t1, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(115, 45, "1", cex = 1)
text(32, 90, "2", cex = 1)
text(-105, 25, "3", cex = 1)
text(-135, -70, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t1 - scale, cex = 0.6)
text(mean(scale), -40, "time (mya)", cex = 0.7)
nodelabels(bg = "white", cex = 0.5)
## Standard tree
plot(ASER[[sample(i, 1)]], colors = cols, ftype = "off")
## Standard tree
plot(ASER[[sample(i, 1)]], colors = cols, ftype = "off")
nodelabels(bg = "white", cex = 0.5)
## Ancestral character state reconstruction without mixed foraging species
colors2 <- c("skyblue", "purple")
cols2 <- setNames(colors2[1:length(unique(forg2))], sort(unique(forg2)))
colorbars2 <-setNames(c("skyblue", "purple"), c("sit and wait", "widely foraging"))
ss2 <- getStates(ASSYM2, "tips")
barcols2 <- setNames(sapply(ss2, function(x, y) y[which(names(y) == x)], y = colorbars2), names(ss2))
i2 <- 1:1000
plotTree.wBars(ASSYM2[[sample(i2, 1)]], log1p(anc_st_rp2_bars), type = "fan", method = "plotSimmap", colors = cols2, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols2, border = FALSE)
t3 <- max(nodeHeights(ctree_rp2))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t3, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(115, 45, "1", cex = 1)
text(45, 85, "2", cex = 1)
text(-105, 35, "3", cex = 1)
text(-137, -68, "4", cex = 1)
text(-45, -65, "5", cex = 1)
text(-8, -95, "6", cex = 1)
text(scale, rep(-23, length(scale)), t3 - scale, cex = 0.9)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd2$ace, piecol = cols2, cex = 0.26)
add.simmap.legend(colors = cols2, prompt = FALSE, x = -120, y = 140, fsize = 0.8)
## Excluding the tuatara
i2T <- 1:1000
colors2T <- c("skyblue", "purple")
cols2T <- setNames(colors2T[1:length(unique(forg2T))], sort(unique(forg2T)))
colorbars2T <-setNames(c("skyblue", "purple"), c("sit and wait", "widely foraging"))
ss2T <- getStates(ASSYMT, "tips")
barcols2T <- setNames(sapply(ss2T, function(x, y) y[which(names(y) == x)], y = colorbars2T), names(ss2T))
plotTree.wBars(ASSYMT[[sample(i2T, 1)]], log1p(anc_st_rp2_bars[-425]), type = "fan", method = "plotSimmap", colors = cols2T, scale = 4, fsize = 0.2, tip.labels = FALSE, part = 0.98, col = barcols2T, border = FALSE)
t4 <- max(nodeHeights(ctree_st2))
tick.spacing <- 25
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(t4, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
text(55, 20, "1", cex = 1)
text(14, 20, "2", cex = 1)
text(-30, 8, "3", cex = 1)
text(-72, -35, "4", cex = 1)
text(-11, -24, "5", cex = 1)
text(-8, -40, "6", cex = 1)
text(scale, rep(-20, length(scale)), t4 - scale, cex = 0.7)
text(mean(scale), -40, "time (mya)", cex = 0.9)
nodelabels(pie = pd2T$ace, piecol = cols2T, cex = 0.26)
add.simmap.legend(colors = cols2T, prompt = FALSE, x = -70, y = 80, fsize = 0.8)